program example_PML
    use common_functions
    use myfem
    use MySparseOperator
    implicit none
    integer, parameter :: sms = 6
    integer, parameter :: nk = 3
    real(kind=8) :: eps = 1.d-10
    real(kind=8) :: O(2), width, height, int_domain(4)
    integer :: m, n, ms_pml(sms), ims, m_pml, n_pml
    integer :: dofs, ii, m1, m2, n1, n2, p
    double precision :: ks(nk)

    real(kind=8) :: k, sigma1, sigma2
    real(kind=8) :: errors(4)
    complex(kind=8), dimension(:), allocatable :: u
    integer :: err

    p = 0

    O = (/ -0.24d+1, -0.22d+1 /)
    width = 4.8d+0
    height = 4.4d+0

    int_domain = (/ -2.d0, 2.d0, -2.d0, 2.d0 /)


    ks = 2*pi* (/ 6, 25, 50 /)
    ms_pml = (/ 8, 16, 32, 64, 128, 256 /)

    write(*,*) 'O', O
    write(*,*) 'width', width, '; height', height
    write(*,*) 'sigma1', sigma1, '; sigma2', sigma2

    do ii = 1, nk, 1
        k = ks(ii)
        sigma1 = 5.d+1 / k
        sigma2 = 5.d+1 / k
        do ims = 1, sms, 1
            m_pml = ms_pml(ims)
            n_pml = m_pml

            m = 12*m_pml
            n = 11*n_pml

            dofs = (m+1)*(n+1)
            allocate(u(dofs))
            call mf_linear_cartesian_DB_zsolver(O, width, height, m, n, a_fun, c_fun, f_fun, g_fun, u)

            m1 = m_pml + 1
            m2 = m - m_pml - 1
            n1 = n_pml + 1
            n2 = n - n_pml - 1
            errors = (/ 1.d0,   &   ! get the L2 norm of u_fun
                &       1.d0,   &   ! get the L2 error between u and  u_fun
                &       1.d0,   &   ! get the H1 semi norm of u_fun
                &       1.d0 /)     ! get the H1-semi errors bet
            call mf_linear_cartesian_zerror(O, width, height, m, n, u_funv, u, errors, m1, m2, n1, n2)
            write(*,'(f10.2,2I5,5f14.7)') k, m, n, errors, errors(4)/errors(3)

            if (allocated(u)) deallocate(u)
        end do
    end do

    if (allocated(u)) deallocate(u)
contains

    subroutine a_fun(x,y,a)
        real(kind=8), intent(in) :: x(:), y(:)
        complex(kind=8), intent(inout) :: a(:,:)

        call a_fun_full(x,y,int_domain,a)
    end subroutine a_fun

    subroutine c_fun(x,y,c)
        real(kind=8), intent(in) :: x(:), y(:)
        complex(kind=8), intent(inout) :: c(:)
        call c_fun_full(x,y,int_domain,c)
    end subroutine c_fun


    subroutine f_fun(x,y,f)
        real(kind=8), intent(in) :: x(:), y(:)
        complex(kind=8), intent(inout) :: f(:)

        real(kind=8), dimension(:), allocatable :: r
        integer :: err, sx, ii
        sx = size(x)
        allocate(r(sx))
        r = dsqrt(x**2+y**2)
        do ii = 1, sx, 1
            if ( r(ii).le.eps ) then
                f(ii) = 0.d0
            elseif ( r(ii).lt.1.d0 .and. r(ii).gt.eps ) then
                f(ii) = -((((36.d0*r(ii)+75.d0)*r(ii)-192.d0)*r(ii)+81.d0)*r(ii))*dcmplx(bessel_j0(k*r(ii)),bessel_y0(k*r(ii)))
                f(ii) = f(ii) + (((((12.d0*k)*r(ii)+(30.d0*k))*r(ii)-(96.d0*k))*r(ii)+(54.d0*k))*r(ii)**2)*dcmplx(bessel_jn(1,k*r(ii)),bessel_yn(1,k*r(ii)))
            else
                f(ii) = 0.d0
            end if
        end do
        if (allocated(r)) deallocate(r)
    end subroutine f_fun

    subroutine g_fun(x,y,v)
        double precision, intent(in) :: x(:), y(:)
        complex(kind=8), intent(inout) :: v(:)
        v = 0.d0
    end subroutine g_fun












    subroutine u_funv(x,y,u,pxy)
        real(kind=8), intent(in) :: x(:), y(:)
        complex(kind=8), intent(inout) :: u(:)
        character, intent(in), optional :: pxy

        real(kind=8), allocatable :: r(:)
        integer :: err,sx
        sx = size(x)
        allocate(r(sx))
        r = dsqrt(x**2+y**2)
        if ( .not.present(pxy) ) then
            where( r.le.eps )
            u = 0.d0
            end where
            where( r.le.1.d0 .and. r.gt.eps )
                u = r**3*(r*(r*(r+3.d0)-12.d0)+9.d0) * &
                    &   dcmplx(bessel_j0(k*r), bessel_y0(k*r))
            end where
            where(r.gt.1.d0)
                u = dcmplx(bessel_j0(k*r), bessel_y0(k*r))
            end where
            return
        end if
        if ( pxy.eq.'x' ) then
            where ( r.lt.eps )
                u = 0.d0
            end where
            where ( r.lt.1.d0 .and. r.ge.eps )
                u = ( r**2 * ( 0.27d+2 + r*( -0.48d+2 + r*( 0.15d+2 + 0.6d+1*r)))*dcmplx(bessel_j0(k*r),bessel_y0(k*r))    &
                    &   - k*r**3 * (0.9d+1 + r*( -0.12d+2 + r*( 0.3d+1 +r)))*dcmplx(bessel_jn(1,k*r),bessel_yn(1,k*r)) )    &
                    &   * (-x/r)
            end where
            where ( r.ge.1.d0  )
                u = k * dcmplx(bessel_jn(1,k*r),bessel_yn(1,k*r)) * x/r
            end where
            u = - u
        end if
        if ( pxy.eq.'y' ) then
            where ( r.lt.eps )
                u = 0.d0
            end where
            where ( r.lt.1.d0 .and. r.ge.eps  )
                u = ( r**2 * ( 0.27d+2 + r*( -0.48d+2 + r*( 0.15d+2 + 0.6d+1*r)))*dcmplx(bessel_j0(k*r),bessel_y0(k*r))    &
                    &   - k*r**3 * (0.9d+1 + r*( -0.12d+2 + r*( 0.3d+1 +r)))*dcmplx(bessel_jn(1,k*r),bessel_yn(1,k*r)) )    &
                    &   * (-y/r)
            end where
            where ( r.ge.1.d0  )
                u = k * dcmplx(bessel_jn(1,k*r),bessel_yn(1,k*r)) * y/r
            end where
            u = - u
        end if

        if (allocated(r)) deallocate(r)
    end subroutine u_funv


    subroutine a_fun_full(x,y,int_domain,a)
        real(kind=8), intent(in) :: x(:), y(:)
        real(kind=8), intent(in) :: int_domain(4)
        complex(kind=8), intent(inout) :: a(:,:)

        integer :: ii, jj, sx

        sx = size(x)

        do ii= 1, sx, 1
            if ( x(ii).gt.int_domain(1) .and. x(ii).lt.int_domain(2) .and. y(ii).gt.int_domain(3) .and. y(ii).lt.int_domain(4) ) then
                a(1,ii) = 1.d0
                a(2,ii) = 1.d0
            elseif ( x(ii).le.int_domain(1) .and. y(ii).gt.int_domain(3) .and. y(ii).lt.int_domain(4) ) then
                a(2,ii) = dcmplx(1.d0,sigma1)
                a(1,ii) = 1.d0/a(2,ii)
            elseif ( x(ii).ge.int_domain(2) .and. y(ii).gt.int_domain(3) .and. y(ii).lt.int_domain(4) ) then
                a(2,ii) = dcmplx(1.d0,sigma1)
                a(1,ii) = 1.d0/a(2,ii)
            elseif ( y(ii).ge.int_domain(4) .and. x(ii).gt.int_domain(1) .and. x(ii).lt.int_domain(2) ) then
                a(1,ii) = dcmplx(1.d0,sigma2)
                a(2,ii) = 1.d0/a(1,ii)
            elseif ( y(ii).le.int_domain(3) .and. x(ii).gt.int_domain(1) .and. x(ii).lt.int_domain(2) ) then
                a(1,ii) = dcmplx(1.d0,sigma2)
                a(2,ii) = 1.d0/a(1,ii)
            elseif ( x(ii).le.int_domain(1) .and. y(ii).le.int_domain(3) ) then
                a(1,ii) = dcmplx(1.d0,sigma2) / dcmplx(1.d0,sigma1)
                a(2,ii) = 1.d0/a(1,ii)
            elseif ( x(ii).ge.int_domain(2) .and. y(ii).le.int_domain(3) ) then
                a(1,ii) = dcmplx(1.d0,sigma2) / dcmplx(1.d0,sigma1)
                a(2,ii) = 1.d0/a(1,ii)
            elseif ( x(ii).ge.int_domain(2) .and. y(ii).ge.int_domain(4) ) then
                a(1,ii) = dcmplx(1.d0,sigma2) / dcmplx(1.d0,sigma1)
                a(2,ii) = 1.d0/a(1,ii)
            elseif ( x(ii).le.int_domain(1) .and. y(ii).ge.int_domain(4) ) then
                a(1,ii) = dcmplx(1.d0,sigma2) / dcmplx(1.d0,sigma1)
                a(2,ii) = 1.d0/a(1,ii)
            end if
        end do
    end subroutine a_fun_full

    subroutine c_fun_full(x,y,int_domain,c)
        real(kind=8), intent(in) :: x(:), y(:)
        real(kind=8), intent(in) :: int_domain(4)
        complex(kind=8), intent(inout) :: c(:)

        integer :: ii, jj, sx

        sx = size(x)

        do ii= 1, sx, 1
            if ( x(ii).gt.int_domain(1) .and. x(ii).lt.int_domain(2) .and. y(ii).gt.int_domain(3) .and. y(ii).lt.int_domain(4) ) then
                c(ii) = - k**2
            elseif ( x(ii).le.int_domain(1) .and. y(ii).gt.int_domain(3) .and. y(ii).lt.int_domain(4) ) then
                c(ii) = (-k**2)*dcmplx(1.d0,sigma1)
            elseif ( x(ii).ge.int_domain(2) .and. y(ii).gt.int_domain(3) .and. y(ii).lt.int_domain(4) ) then
                c(ii) = (-k**2)*dcmplx(1.d0,sigma1)
            elseif ( y(ii).ge.int_domain(4) .and. x(ii).gt.int_domain(1) .and. x(ii).lt.int_domain(2) ) then
                c(ii) = (-k**2)*dcmplx(1.d0,sigma2)
            elseif ( y(ii).le.int_domain(3) .and. x(ii).gt.int_domain(1) .and. x(ii).lt.int_domain(2) ) then
                c(ii) = (-k**2)*dcmplx(1.d0,sigma2)
            elseif ( x(ii).le.int_domain(1) .and. y(ii).le.int_domain(3) ) then
                c(ii) = (-k**2)*dcmplx(1.d0,sigma2) * dcmplx(1.d0,sigma1)
            elseif ( x(ii).ge.int_domain(2) .and. y(ii).le.int_domain(3) ) then
                c(ii) = (-k**2)*dcmplx(1.d0,sigma2) * dcmplx(1.d0,sigma1)
            elseif ( x(ii).ge.int_domain(2) .and. y(ii).ge.int_domain(4) ) then
                c(ii) = (-k**2)*dcmplx(1.d0,sigma2) * dcmplx(1.d0,sigma1)
            elseif ( x(ii).le.int_domain(1) .and. y(ii).ge.int_domain(4) ) then
                c(ii) = (-k**2)*dcmplx(1.d0,sigma2) * dcmplx(1.d0,sigma1)
            end if
        end do
    end subroutine c_fun_full









end program example_PML